/*******************************************
  * Producted by Earthson		           *
  * Email: Earthson.Landy@gmail.com        *
  * http://www.earthson.net		           *
  ******************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <time.h>
#include "numbertheory.h"
/**
  * The default type of number is int, you may need
  * other types(eg. long long) to deal with much more larger number^^
  **/


const int N = 1000000;
int ph[N+1] = {0};     	///phi(x) for x<=N
int prime[N+1];  		///prime table,prime[0] is the number of primes in table^^
int sig[N+1] = {0};     ///sigma(x) for x<=N,sigma(x) is sum of x's all productors^^
int productorcnt[N+1]; 	///count of productor numbers frome 1 ro N

/**
  * Eculid Algorithm-----GCD
  **/
int gcd (int a, int b)
{
    return b ? gcd (b, a % b) : a;
}

/**
  * LCM of k numbers
  **/
int lcm (int *table, int k)
{
    int result = 1;
    for (int i = 0; i < k; i++)
        result *= table[i] / gcd (result, table[i]); ///may overfow...
    return result;
}

/**
  * Extended Euclid's Algorithm
  * ax+by=gcd(a,b);
  **/
int extended_gcd (int a, int b, int &x, int &y)
{
    if (!b)
        return x = 1, y = 0, a;
    int ret = extended_gcd (b, a % b, x, y), tmp = x;
    x = y, y = tmp - (a / b) * y;
    return ret;
}

/**
  * a*b mod m
  * you may need this to deal with big a&b that a*b overflows, this function helps you do that^^
  * m*2<INT_MAX
  **/
int a_Multiply_b_Mode_m (int a, int b, int m)
{
    int result = 0;
    b %= m;
    while (a %= m, b)
    {
        if (b & 1)
            result += a, result %= m;
        a <<= 1, b >>= 1;
    }
    return result;
}
/**
  * Successive square method
  * m*m < INT_MAX
  * O(lgn)
  **/
int Power_k_Mod_m (int x, int k, int m) ///x^k(mod m)
{
    int result = 1;
    while (x %= m, k)
    {
        if (k & 1)
            result *= x, result %= m;
        x *= x, k >>= 1;
    }
    return result;
}
int power (int x, int k)
{
    int result = 1;
    while (k)
    {
        if (k & 1)
            result *= x;
        x *= x, k >>= 1;
    }
    return result;
}

/**
  * Computing k-th Roots Modulo m
  * x^k = b (mod m) ,return x
  * Be sure that gcd(b,m)=1,and gcd(k,phi(m))=1.
  **/
int Root_k_Mode_m (int b, int k, int m)
{
    int u, v, phm;
    phm = phi (m); ///need a function to get phi(m)
    v = extended_gcd (k, m, u, v);
    if (u < 0)
        u += phm / v;
    return Power_k_Mod_m (b, u, m);
}

/**
  * 1+q+q^2+...+q^n-1
  * O(logn)
  **/
int sum_Mode_m (int q, int n, int m) ///(1+q+q^2+...+q^n-1) mod m
{
    int result = 0, temp = 1;
    while (q %= m, n)
    {
        if (n % 2)
            result *= q, result += temp, result %= m;
        temp *= (1 + q), q *= q, temp %= m, n >>= 1;
    }
    return result;
}

/**
  * Prime number sieve Method
  * O(n)
  **/
int get_prime() ///get primes from 1 to N in liner time
{
    memset (prime, 0, sizeof (int) * (N + 1));
    for (int i = 2; i <= N; i++)
    {
        if (!prime[i])
            prime[++prime[0]] = i;
        for (int j = 1; j <= prime[0] && prime[j] <= N / i; j++)
        {
            prime[prime[j]*i] = 1;
            if (i % prime[j] == 0) break;
        }
    }
    return prime[0];
}

/**
  * Miller-Rabin primality test
  **/
bool MR_primality_test (int x)	///x*x<INT_MAX
{
    bool flag = true;
    const int cnt = 50;
    srand ( (unsigned) time (NULL));
    int a, k = 0, q = x - 1, s;
    while (q % 2 == 0) k++, q /= 2;
    for (int i = 0; i < cnt && flag; i++)
    {
        a = rand() % (x - 1), a++;
        s = Power_k_Mod_m (a, q, x); ///need Power_k_Mod_m
        if (s == 1) continue;
        for (int j = k; j && flag; j--)
        {
            if (s == x - 1)
                flag = false;
            s *= s, s %= x;
        }
        flag = flag ? false : true;
    }
    return flag;
}

bool isPrime (int x)
{
    return MR_primality_test (x);
}
/**
  * phi(n) !!!No prime table
  **/
int phi (int n)
{
    int result = n;
    for (int i = 2; i <= n / i; i++)
    {
        if (n % i == 0)  ///i is always a prime^^
            result /= i, result *= i - 1;
        while (n % i == 0)
            n /= i;
    }
    if (n > 1)
        result /= n, result *= n - 1;
    return result;
}

/**
  * phi(n) !!!need prime table
  * if size of table is N, be sure that n*n<=N
  **/
int phi_2 (int n)
{
    int result = n;
    for (int i = 1; i <= prime[0] && prime[i] <= n; i++)
    {
        if (n % prime[i] == 0)
            result /= prime[i], result *= prime[i] - 1;
        while (n % prime[i] == 0)
            n /= prime[i];
    }
    if (n > 1)
        result /= n, result *= n - 1;
    return result;
}

/**
  * phi from 1 to N
  * you can also get a prime table frome this function
  * O(n)
  **/
void phi_all() ///get phi(n) from 1 to N in liner time
{
    memset (ph, 0, sizeof (int) * (N + 1));
    prime[0] = 0;
    for (int i = 2; i <= N; i++)
    {
        if (!ph[i])
            prime[++prime[0]] = i, ph[i] = i - 1;
        for (int j = 1; prime[j] <= N / i && j <= prime[0]; j++)
        {
            if (i % prime[j])
                ph[prime[j]*i] = ph[prime[j]] * ph[i];
            else
            {
                ph[prime[j]*i] = ph[i] * prime[j];
                break;
            }
        }
    }
}

/**
  * sigma(n) ||| need prime table!
  * sum of all productors of n
  **/
int sigma (int n)
{
    int pk = 1, result = 0;
    for (int i = 1; prime[i] <= n / prime[i] && i <= prime[0]; i++)
    {
        while (n % prime[i] == 0)
            n /= prime[i], pk *= prime[i], pk++;
        result *= pk, pk = 1;
    }
    if (n > 1)
        pk *= n, pk++, result *= pk;
    return result;
}

/**
  * sigma(n) ||| if you do not have a prime table^^
  * sum of all productors of n
  **/
int sigma2 (int n)
{
    int pk = 1, result = 0;
    for (int i = 2; i <= n / i; i++)
    {
        while (n % prime[i] == 0) n /= i, pk *= i, pk++;
        result *= pk, pk = 1;
    }
    if (n > 1)
        pk *= n, pk++, result *= pk;
    return result;
}

/**
  * sigma(n) ||| also can get a prime table^^
  * sum of all productors of n for all n no more than N^^
  * O(n)
  **/
int min_pk[N+1];
void sigma_all()
{
    memset (sig, 0, sizeof (int) * (N + 1));
    sig[1] = 1, min_pk[1] = 1, prime[0] = 0;
    for (int i = 2; i <= N; i++)
    {
        if (!sig[i])
            prime[++prime[0]] = i, sig[i] = i + 1, min_pk[i] = i + 1;
        for (int j = 1; prime[j] <= N / i && j <= prime[0]; j++)
        {
            if (i % prime[j])
            {
                sig[i*prime[j]] = sig[i] * sig[prime[j]];
                min_pk[i*prime[j]] = min_pk[prime[j]];
            }
            else
            {
                min_pk[i*prime[j]] = min_pk[i] * prime[j] + 1;
                sig[i*prime[j]] = (sig[i] / min_pk[i]) * min_pk[i*prime[j]];
                break;
            }
        }
    }
}

/**
  * productors count of all number from 1 to N
  * it can get a prime table too
  * O(n)
  **/
int min_cnt[N+1];
void productor_all()
{
    memset (productorcnt, 0, sizeof (int) * (N + 1));
    productorcnt[1] = 1, min_cnt[1] = 1, prime[0] = 0;
    for (int i = 2; i <= N; i++)
    {
        if (!productorcnt[i])
            prime[++prime[0]] = i, productorcnt[i] = 2, min_cnt[i] = 2;
        for (int j = 1; prime[j] <= N / i && j <= prime[0]; j++)
        {
            if (i % prime[j])
            {
                productorcnt[i*prime[j]] = productorcnt[i] * productorcnt[prime[j]];
                min_cnt[i*prime[j]] = min_cnt[prime[j]];
            }
            else
            {
                min_cnt[i*prime[j]] = min_cnt[i] + 1;
                productorcnt[i*prime[j]] = (productorcnt[i] / min_cnt[i]) * min_cnt[i*prime[j]];
                break;
            }
        }
    }
}
/**
  * get n! mod p,p is prime number
  * almost O(n/ln(n))
  * You need a prime table^^
  * if n! = a * p^k,and a mod p != 0
  * then return a mod m,and you can get k from the Variable:over^^
  **/
int factorial_mod_p (int n, int p, int &over) ///p <= prime[prime[0]]
{
    if(p == 1) return 0;
    int result = 1, cnt = 0, temp;
    for(int i = 1; i < prime[0] && prime[i] <= n; i++)
    {
        temp = n, cnt = 0;
        while(temp)
            temp /= prime[i], cnt += temp;
        temp = prime[i];
        if(prime[i] == p)
            over = cnt, cnt = 0;
        while(cnt)
        {
            if(cnt & 1)
                result *= temp, result %= p;
            temp *= temp, cnt >>= 1, temp %= p;
        }
    }
    return result;
}
///O(n)
int factorial_mod_p_liner(int n, int p, int &over)
{
    if(p == 1) return 0;
    int result = 1;
    over = 0;
    for(int i = 2; i <= n; i++)
    {
        result *= i;
        while(result % p == 0)
            result /= p, over++;
        result %= p;
    }
    return result;
}

/**
  * Combinations( C(a, b) ) mod m
  * O(b^2*lgb) m*m<INT_MAX
  **/
int combination_Mod_m (int a, int b, int m)
{
    const int gate = 2000000000;  ///larger gate for more Optimization
    if (b > a) return 0;
    b = (a - b < b) ? a - b : b;
    int d[b], tmp = 1 ,l = 0, h, j, g;
    for (l = 0, j = a - b + 1; a - l >= j; l++)
    {
        d[l] = a - l;
        while(gate / d[l] >= j && a - l != j)
            d[l] *= j, j++;
    }
    for (j = 2, h = b; j <= h; h--)
    {
        tmp *= h;
        while (gate / tmp >= j && h != j)
            tmp *= j, j++;
        for (int i = 0; tmp != 1; i++)
        {
            g = gcd (tmp, d[i]);
            d[i] /= g, tmp /= g;
        }
    }
    int result = 1;
    for (int i = 0; i < l; i++)
        d[i] %= m, result *= d[i], result %= m;
    return result;
}

///O(m)
///p is prime number, p*p < INT_MAX
int combination_Mod_p (int n, int m, int p)
{
    if (m > n) return 0;
    m = (n - m < m) ? n - m : m;
    int a = 1, b = 1, x, y, pcnt = 0;
    for (int i = 1; i <= m; i++)
    {
        a *= n - i + 1, b *= i;
        while (a % p == 0) a /= p, pcnt++;
        while (b % p == 0) b /= p, pcnt--;
        b %= p, a %= p;
    }
    if (pcnt) return 0;
    extended_gcd(b, p, x, y);
    if (x < 0) x += p;
    x *= a, x %= p;
    return x;
}

///O(m)
///p is prime and p^k < INT_MAX
int combination_Mod_pk (int n, int m, int p, int k)
{
    if (m > n) return 0;
    m = (n - m < m) ? n - m : m;
    int a = 1, b = 1, x, y, pa = 1, pb = 1, pcnt = 0;
    int q = power(p, k);
    for (int i = 1; i <= m; i++)
    {
        a *= n - i + 1, b *= i;
        while(a % p == 0) pa *= p, a /= p;
        while(b % p == 0) pb *= p, b /= p;
        while (pa % q == 0) pa /= q, pcnt++;
        while (pb % q == 0) pb /= q, pcnt--;
        b %= q, a %= q;
    }
    if(pa < pb) pcnt--, pa *= q;
    pa /= pb, a *= pa, a %= q;
    if (pcnt) return 0;
    extended_gcd(b, q, x, y);
    if (x < 0) x += q;
    x *= a, x %= q;
    return x;
}
/// s*s < INT_MAX , and s can be any positive number
/// O(m*lg(n))
int combination_Mod_s (int n, int m, int s)
{
    if (m > n) return 0;
    m = (n - m < m) ? n - m : m;
    int a = 1, b = 1, x, y, g, pcnt = 0;
    int pa = 1, pb = 1; ///To be sure that pa and pb is not overflow
    for (int i = 1; i <= m; i++)
    {
        a *= n - i + 1, b *= i;
        while( (g = gcd(a, s)) > 1) pa *= g, a /= g;
        while( (g = gcd(b, s)) > 1) pb *= g, b /= g;
        g = gcd(pa, pb), pa /= g, pb /= g;
        while (pa % s == 0) pa /= s, pcnt++;
        while (pb % s == 0) pb /= s, pcnt--;
        b %= s, a %= s;
    }
    while(pcnt) pcnt--, pa *= s;
    pa /= pb, pa %= s, a *= pa, a %= s;
    if (pcnt) return 0;
    extended_gcd(b, s, x, y);
    if (x < 0) x += s;
    x *= a, x %= s;
    return x;
}
///need prime table,and n < N
///almost less than O(n)
int combination_Mod_h (int n, int m, int h)
{
    if(h == 1) return 0;
    int result = 1, cnt = 0, temp;
    for(int i = 1; i < prime[0] && prime[i] <= n; i++)
    {
        temp = n, cnt = 0;
        while(temp)
            temp /= prime[i], cnt += temp;
        temp = n - m;
        while(temp)
            temp /= prime[i], cnt -= temp;
        temp = m;
        while(temp)
            temp /= prime[i], cnt -= temp;
        temp = prime[i];
        while(cnt)
        {
            if(cnt & 1)
                result *= temp, result %= h;
            temp *= temp, cnt >>= 1, temp %= h;
        }
        if(result == 0) return 0;
    }
    return result;
}

/**
  * Lucas' theorem---combination number mod prime p
  * Be sure that p*p<INT_MAX
  **/
int lucas (int m, int n, int p)
{
    int result = 1;
    while (m && n && result)
    {
        result *= combination_Mod_p (m % p, n % p, p);
        result %= p, m /= p, n /= p;
    }
    return result;
}

/**
  * small Combinations
  * C(n,m)<INT_MAX
  * almost O(m)
  **/
int combination (int n, int m)
{
    int result = 1, g;
    if (m > n) return 0;
    m = (n - m < m) ? n - m : m;
    for (int i = 0; i < m; i++)
    {
        if (result % (i + 1) == 0)
            result /= i + 1, result *= n - i;
        else if ( (n - i) % (i + 1) == 0)
            result *= (n - i) / (i + 1);
        else
        {
            g = gcd (result, i + 1);
            result /= g;
            result *= (n - i) / ( (i + 1) / g);
        }
    }
    return result;
}

/**
  * Linear congruence theorem
  * x=b(mod m)
  * x=c(mod n)
  * gcd(m,n)=1,,0<=x<mn
  **/
int linear_congruence_theorem (int b, int c, int m, int n)
{
    int x, y;
    y = extended_gcd (m, n, x, y);
    while (c < b) c += n / y;
    x *= c - b, x = m * x + b, x %= m * n;
    return x;
}


/**
  * Polya for Dihedral group(n rotational symmetries and n reflection symmetries)
  * n points and c colors
  **/
int polya (int c, int n)
{
    int result = 0;
    for (int i = 0; i < n; i++) result += power (c, gcd (i, n));
    if (n & 1)
        for (int i = 0; i < n; i++)
            result += power (c, n / 2 + 1);
    else
        for (int i = 0; i < n / 2; i++)
            result += power (c, n / 2) + power (c, n / 2 + 1);
    return result;
}


/**
  * Baby step, Giant step for Discrete Logarithm
  * a^x = b (mod p), get the value of x
  * O(sqrt(n))
  **/
int giantStep, tableSize;
int *hashTable, *toHashTable, *linkTable;
void BSGS_init(int p)
{
	if(hashTable != NULL)
		delete [] hashTable;
	if(toHashTable != NULL)
		delete [] toHashTable;
	if(linkTable != NULL)
		delete [] linkTable;
	tableSize = giantStep = ceil(sqrt(p));
	while(!isPrime(tableSize)) tableSize++;
	hashTable = new int[tableSize];
	toHashTable = new int[giantStep];
	linkTable = new int[giantStep];
	for(int i = 0; i < tableSize; i++)
		hashTable[i] = -1;
}

int hash(int tohash)
{
	return tohash % tableSize;
}

int hashFind(int tohash)
{
	int pos = hash(tohash);
	int linker = hashTable[pos];
	while(linker != -1)
	{
		if(toHashTable[linker] == tohash)
			return linker;
		linker = linkTable[linker];
	}
	return -1;
}

void getHashTable(int a, int p)
{
	int tmp = 1, pos;
	a %= p;
	for(int i = 0; i < giantStep; i++)
	{
		toHashTable[i] = tmp;
		pos = hash(tmp);
		if(hashFind(tmp) == -1)
		{
			linkTable[i] = hashTable[pos];
			hashTable[pos] = i;
		}
		tmp *= a, tmp %= p;
	}
}

int babyStep_gaintStep(int a, int b, int p)
{
	BSGS_init(p);
	getHashTable(a, p);
	int step = Power_k_Mod_m(a, p - giantStep - 1, p);
	int tmp = b, ret;
	tmp %= p;
	for(int i = 0; i < giantStep; i++)
	{
		ret = hashFind(tmp);
		if( ret != -1)
			return i * giantStep + ret;
		tmp *= step, tmp %= p;
	}
	return -1;
}

/**
  * Continued Fraction
  **/
int fractiontable[30] = { 0 };
void continuedFraction(int a, int b)
{
	if(b)
	{
		fractiontable[++fractiontable[0]] = a/b;
		continuedFraction(b, a%b);
	}
}
